
#like 1.1, but use uneven weight matrix and negative values of \rho

rm(list=ls())
library(epicalc)
zap()
#Spatial spec correct, omitted variable in X\betas (correlated with \rho)

sar.f <- function(es, x, y, w){
    rho <- es[1]
    s2 <- es[2]^2
    if (s2<0) return(-1e20)
    k <- ncol(x)
    b <- as.matrix(es[3:(k+2)])
    llik <- 0
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    y <- as.matrix(y)
    
    for(i in 1:d){
        A <- diag(1, nd, nd)-rho*w[,,i]
        dA <- det(A)
        if (dA <=0){
            llik <-  -1e20 
            break
        } else {
            aux <- log(dA) -nd/2 * log(s2)- (1/(2*s2)) * t(A %*% y[((i-1)*nd+1):((i-1)*nd+nd),] - as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b) %*% (A %*% y[((i-1)*nd+1):((i-1)*nd+nd),]- as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b)
            llik <- llik + aux
        }
    }
    return(llik)
}

sar <- function(x, y, w, stval,vars=c(""), routit=100000){
    mod <- optim(stval, sar.f, hessian = TRUE, method = "BFGS", 
            control = list(fnscale = -1, REPORT = 10, 
            reltol = 1e-16, trace = 1, 
            maxit = routit), x=x, y=y, w=w)
    vacov <- as.matrix(solve(-1*mod$hessian, tol=1e-24))
    se <- as.matrix(sqrt(diag(vacov)))
    zstat <- as.matrix(mod$par / se)
    pvalues <- round(2*pnorm(abs(zstat),lower.tail=FALSE),4)
    
    cat("results in order: rho, s, b \n")
    cat("parameter estimates: ", mod$par, "\n")
    cat("standard error: ", se, "\n")
    cat("z-statistics: ", zstat, "\n")
    cat("p-values: ", pvalues, "\n")
    list(llik=mod$value, par=mod$par, se=se, zstat=zstat, 
        pvalues = pvalues, 
        hessian=mod$hessian, vacov=vacov
        )
}



genY <- function(nd, d, sig2, rho, b0, b1, b2, w.cor){ #Generate data from correct DGP

               
    ytrue <- rep(NA, d*nd)
    xother1 <- rep(NA, d*nd)        
    xother2 <- rep(NA, d*nd)      
    xotherpanel <- rnorm(d)
      

    for (j in 1:d){
        iA <- solve(diag(1,nd)-rho * w.cor[,,j])
        xotheraux1 <- rnorm(nd)
        xotheraux2 <- rep(xotherpanel[j], nd)
        yaux <- iA %*% (b0 + b1 * xotheraux1 + b2 * xotheraux2) + iA%*%rnorm(nd, sd=sqrt(sig2))
        ytrue[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux
        xother1[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux1
        xother2[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux2
    }
    return(list(ytrue=ytrue, xother1=xother1, xother2=xother2, w.cor=w.cor))
}


#repetitions
r <- 1000
    
#specify parameters
sig2 <- 1.0
rho <- -.8
b0 <- .5
b1 <- 1
b2 <- 3
sca <- 1

nd<-10
d<-100
#set.seed(021175)
set.seed(08202002)

#initiate correct connectivity matrix
w.rd <- array(NA, c(nd, nd, d))
row.st <- function(rn, w){
    w[rn,]/sum(w[rn,])
}
for (i in 1:d){
    wblock <- matrix(runif(nd*nd), nd, nd)
    diag(wblock) <- 0
    wblock <- t(sapply(1:nd, row.st, w=wblock))
    w.rd[,,i] <- wblock
}    
   
#initate naive even connectivit matrix

wall.ev <- array(NA, c(nd, nd, d))
wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
diag(wraw) <- 0
wall.ev[,,1:d] <- wraw


mod.true <- vector("list", r)
mod.lack <- vector("list", r)
i <- 1


while (i <= r){    
    mcdat <- genY(nd, d, sig2, rho, b0, b1, b2, w.rd)
    xfull <- cbind(1, mcdat$xother1, mcdat$xother2)
    xlack <- cbind(1, mcdat$xother1)
    y <- mcdat$ytrue
    stvfull <- runif((ncol(xfull) + 2))*.1
    stvlack <- runif((ncol(xlack) + 2))*.1
    mod.true[[i]] <- try(sar(xfull, y*sca, wall.ev, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.true[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvfull <- runif((ncol(xfull) + 2))*.1
            mod.true[[i]] <- try(sar(xfull, y*sca, wall.ev, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF TRUE MODEL\n")
                break
            }
    }    
#    mod.lack[[i]] <- try(sar(xlack, y*sca, wall.ev, stval=stvlack, vars=c("const", "other")), TRUE)
#    ect <- 1 #initiate error counter
#    while(class(mod.lack[[i]]) == "try-error"){
#            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
#            stvlack <- runif((ncol(xlack) + 2))*.1
#            mod.lack[[i]] <- try(sar(xlack, y*sca, wall.ev, stval=stvlack, vars=c("const", "other")), TRUE)
#            ect <- ect+1
#            if (ect == 10) {
#                i <-  r+1
#                cat("\n\n\n\n ROUTINE FAILED BC OF FALSE MODEL\n")
#                break
#            }
 #   }
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

getpars <- function(i, model){
    model[[i]]$par
}
getllik <- function(i, model){
    model[[i]]$llik
}
pars.true <- unlist(sapply(1:1000, getpars, mod.true), recursive=FALSE)
pars.true <- matrix(pars.true, r, length(stvfull), byrow=TRUE)

#pars.lack <- unlist(sapply(1:1000, getpars, mod.lack), recursive=FALSE)
#pars.lack <- matrix(pars.lack, r, length(stvlack), byrow=TRUE)

llik.true <- unlist(sapply(1:r, getllik, mod.true), recursive=FALSE)
#llik.lack <- unlist(sapply(1:r, getllik, mod.lack), recursive=FALSE)

quantile(pars.true[,1], c(.025,.5,.975))
#quantile(pars.lack[,1], c(.025,.5,.975))

bias.true <- mean(pars.true[,1]-rho)
#bias.lack <- mean(pars.lack[,1]-rho)

bias.true
#bias.lack

mean(llik.true)
mean(llik.lack)


#true model has better fit:
#ttest of likelihood
t.test(llik.true,llik.lack,paired=TRUE)

#llratio test of mean llik
dgf=1
r <- -2 * (mean(llik.lack)-mean(llik.true))
p <- 1-pchisq(r, dgf)
cat("degrees of freedom: ", dgf, "\n")
cat("likelihood ratio statistic: ", r,"\n")
cat("p value: ", p, " (under H0 that the restricted model is correct)\n\n")
# number of lliks greater for full model
check<- llik.true>llik.lack
sum(check)
#setwd("/Users/datalabuser/Documents/Project Lead Donor/Analysis/R&R.IO")
setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
#save(mod.lack, file="mcnold.w_wrong.plus3.omit.2.out")
save(mod.true, file="mcnold.w_wrong.RNeg0.8.BetaPlus3.full.out")
#rm(list=ls())


r=1000


##############################################################
###Now re-run analysis with correct betas AND correct W matrix
mod.true <- vector("list", r)
i <- 1

while (i <= r){    
    mcdat <- genY(nd, d, sig2, rho, b0, b1, b2, w.rd)
    xfull <- cbind(1, mcdat$xother1, mcdat$xother2)
    y <- mcdat$ytrue
    stvfull <- runif((ncol(xfull) + 2))*.1
    mod.true[[i]] <- try(sar(xfull, y*sca, w.rd, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.true[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvfull <- runif((ncol(xfull) + 2))*.1
            mod.true[[i]] <- try(sar(xfull, y*sca, w.rd, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF TRUE MODEL\n")
                break
            }
    }    
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

getpars <- function(i, model){
    model[[i]]$par
}
getllik <- function(i, model){
    model[[i]]$llik
}
pars.true <- unlist(sapply(1:1000, getpars, mod.true), recursive=FALSE)
pars.true <- matrix(pars.true, r, length(stvfull), byrow=TRUE)

llik.true <- unlist(sapply(1:r, getllik, mod.true), recursive=FALSE)

quantile(pars.true[,1], c(.025,.5,.975))

bias.true <- mean(pars.true[,1]-rho)

bias.true

mean(llik.true)
#setwd("/Users/datalabuser/Documents/Project Lead Donor/Analysis/R&R.IO")
setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
save(mod.true, file="mcnold.w_true.RNeg0.8.BetaPlus3.full.out")
